Rheological constitutive equation for model of soft glassy materials 
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We solve exactly and describe in detail a simplified scalar model for the low frequency shear 
rheology of foams, emulsions, slurries, etc. [P. SoUich, F. Lequeux, P. Hebraud, M.E. Gates, Phys. 
Rev. Lett. 78, 2020 (1997)]. The model attributes similarities in the rheology of such "soft glassy 
materials" to the shared features of structural disorder and metastability. By focusing on the 
dynamics of mesoscopic elements, it retains a generic character. Interactions are represented by a 
mean-field noise temperature x, with a glass transition occurring at a; = 1 (in appropriate units). 
The exact solution of the model takes the form of a constitutive equation relating stress to strain 
history, from which all rheological properties can be derived. For the linear response, we find that 
both the storage modulus G' and the loss modulus G" vary with frequency as a;^~^ for 1 < a; < 2, 
becoming flat near the glass transition. In the glass phase, aging of the moduli is predicted. The 
steady shear flow curves show power law fluid behavior for x < 2, with a nonzero yield stress in 
the glass phase; the Cox-Merz rule does not hold in this non-Newtonian regime. Single and double 
step strains further probe the nonlinear behavior of the model, which is not well represented by the 
BKZ relation. Finally, we consider measurements of G' and G" at finite strain amplitude 7. Near 
the glass transition, G" exhibits a maximum as 7 is increased in a strain sweep. Its value can be 
strongly overestimated due to nonlinear effects, which can be present even when the stress response 
is very nearly harmonic. The largest strain 7c at which measurements still probe the linear response 
is predicted to be roughly frequency-independent. 

PACS numbers: 83.20.-d, 83.70.Hq, 05.40-l-j. To appear in Physical Review E (July 1998). 



I. INTRODUCTION 

Many soft materials, such as foams, emulsions, pastes 
and slurries, have intriguing rheological properties. Ex- 
perimentally, there is a well-developed phenomenology 
for such systems: their nonlinear flow behavior is often 
fit to the form a = A + 57" where a is shear stress 
and 7 strain rate. This is the Herschel-Bulkeley equa- 
tion HI]; or (for A = Q) the "power-law fluid" ||-|. 
For the same materials, linear or quasi-linear viscoelas- 
tic measurements often reveal storage and loss moduli 
G'{ijj), G"{uo) in nearly constant ratio [G" jG is usually 
about 0.1) with a frequency dependence that is either a 
weak power law (clay slurries, paints, microgels) or neg- 
ligible (tomato paste, dense emulsions, dense multilayer 
vesicles, colloidal glasses) p|-[lO|l. This behavior persists 
down to the lowest accessible frequencies (about 10~^~ 
1 Hz depending on the system), in apparent contradic- 
tion to linear response theory, which requires that G"{uj) 
should be an odd function of lu. This behavior could 
in principle be due to slow relaxation modes below the 
experimentally accessible frequency range (see Fig. ||). 
Each of those would cause a drop in G'{uj) and a bump 
in G" (uj) as the frequency is tracked downward. However, 
where the search for system specific candidates for such 
slow modes has been carried out (for the case of foams 
and dense emulsions, for example, see |1 1| ) , it has not 
yielded viable candidates; it therefore seems worthwhile 
to look for more generic explanations of the observed be- 
havior. 

Indeed, the fact that similar anomalous rheology 
should be seen in such a wide range of soft materials sug- 



gests a common cause. In particular, the frequency de- 
pendence indicated above points strongly to the generic 
presence of slow "glassy" dynamics persisting to arbi- 
trarily small frequencies. This feature is found in several 
other contexts |lj-|lj], such as the dynamics of elastic 
manifolds in random media llsUTq] . The latter is sugges- 
tive of rheology: charge density waves, vortices, contact 
lines, etc. can "flow" in response to an imposed "stress". 

In a previous letter [|l^ it was argued that glassy dy- 
namics is a natural consequence of two properties shared 
by all the soft materials mentioned above: structural dis- 
order and metastability. In such "soft glassy materials" 
(SGMs), thermal motion alone is not enough to achieve 
complete structural relaxation. The system has to cross 
energy barriers (for example those associated with re- 
arrangement of droplets in an emulsion) that are very 
large compared to typical thermal energies. It there- 
fore adopts a disordered, metastable configuration even 
when (as in a monodisperse emulsion or foam) the state 
of least free energy would be ordered (l8|. The impor- 
tance of structural disorder has previously been noted in 
more specific contexts [k 11 19-E3], but its unifying role 
in rheological modeling can be more easily appreciated 
by focusing on the class of SGMs as a whole. 

In Ref. U%, a minimal, scalar model for the generic 
rheology of SGMs was introduced, which incorporates the 
above ideas. We refer to this model as the "soft glassy 
rheology" (SGR) model in the following. The main con- 
tribution of the present publication is the exact solution 
of this model; at the same time, we also provide more 
detailed analytical and numerical support for the results 
announced in IIitI . The exact solution is in the form of a 
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FIG. 1. Sketch of frequency dependence of linear moduli, 
showing possible slow relaxation modes at frequencies below 
the measurement window. 



minimal number of features common to all SGMs, leav- 
ing aside as much system specific detail as possible. One 
important feature is the "glassiness" , i.e., the effects of 
structural disorder and metastability. We model this us- 
ing a fairly intuitive picture of a glass: it consists of local 
"elements" (we will be more specific later about what we 
mean by these in the context of SGMs) which are trapped 
in "cages" formed by their neighbors so that they can- 
not move. Occasionally, however, a rearrangement of the 
elements may be possible, due to thermal activation, for 
example. Glass models of this kind are commonly re- 
ferred to as "trap models" and have been studied by a 
large number of authors (see e.g. Refs. [ p^ , p4| -pO|] ) . An 
alternative to such models would be, for example, mode- 
coupling theories |3iy32|] which, at least in their simplest 
form, neglect all (thermally) activated processes. We pre- 
fer trap models for our purposes, because they are simpler 
and also generally more physically transparent pS] . 



A. Bouchaud's glass model 



constitutive equation relating the (shear) stress at a given 
time to the strain history. We use this to study a range of 
linear and nonlinear rheological properties of the model; 
qualitative comparisons with experimental data show 
that these capture many generic rheological characteris- 
tics of SGMs. We do not attempt more quantitative fits 
to experimental data for specific materials because the 
model in its present form is almost certainly too over- 
simplified for this purpose. We do however hope to carry 
out such a more quantitative study in future work, once 
the remaining ambiguities in the interpretation of the 
model parameters (see Sec. VI) have been cla rified and 
some of the improvements suggested in Sec. VI] have 
been incorporated into the model. 

We introduce the SGR model in Sec. O, along with 
Bouchaud's glass model on which it builds. Sec. Ill con- 



tains our main result, the constitutive equation. Its pre- 
dictions in the linear response regime are discussed in 
Sec. IV, while in Sec. M we analyse several nonlinear sce- 
narios including steady shear flow, shear startup, large 
step strains and large oscillatory strains. The physical 
significance and interpretation of the various parameters 
of the SGR model is not obvious; in Sec. VI we discuss 
in more detail the "noise temperature" x and "attempt 
frequency" Fq of the model. Our results are summarized 
in Sec. IvH. 



II. THE SGR MODEL 

The SGR model is a phenomenological model that aims 
to explain the main features of SGM rheology (both lin- 
ear and nonlinear) as described above. To apply to a 
broad range of materials, such a model needs to be rea- 
sonably generic. It should therefore incorporate only a 



Bouchaud formalized the above intuitive trap picture 
of a glass into a one-element model 1 12 1^ : an individ- 
ual clement "sees" an energy landscape of traps of var- 
ious depths E; when activated, it can "hop" to another 
trap. Bouchaud assumed that such hopping processes are 
due to thermal fluctuations. In SGMs, however, this is 
unlikely as A:bT is very small compared to typical trap 
depths E (see Sec. VI). The SGR model assumes instead 
that the "activation" in SGMs is due to interactions: a 
rearrangement somewhere in the material can propagate 
and cause rearrangements elsewhere. In a mean-field 
spirit, this coupling between elements is represented by 
an effective temperature (or noise level) x. This idea is 
fundamental to the SGR model. 

The equation of motion for the probability of finding 
an element in a trap of depth E at time t is [p^p^p^ 



d_ 
di 



P{E, t) = -Foe-^/- PiE, t) + r(t) p{E) 



(1) 



In the first term on the rhs, which describes elements 
hopping out of their current traps, Fq is an attempt fre- 
quency for hops, and exp(— iJ/x) is the corresponding 
activation factor. The second term represents the state 
of these elements directly after a hop. Bouchaud made 
the simplest possible assumption that the depth of the 
new trap is completely independent of that of the old 
one; it is simply randomly chosen from some "prior" dis- 
tribution of trap depths p{E). The rate of hopping into 
traps of depth E is then p{E) times the overall hopping 
rate, given by 



T{t) = Fo ( e-^/^ ) =ro f dE P{E, t) e"^/^ (2) 



-E/x 



Bouchaud's main insight was that the model (0) can 
describe a glass transition if the density of deep traps 



has an exponential tail, p{E) ~ expf— i?/xg), say. Why 
is this? The steady state of eq. (|l|), if one exists, is 
given by Peq{E) oc exp{E/x)p{E); the Boltzmann fac- 
tor exp(i<^/x) (no minus here because trap depths are 
measured from zero downwards) is proportional to the 
average time spent in a trap of depth E. At x = a;g, 
it just cancels the exponential decay of p{E), and so the 
supposed equilibrium distribution Pcq{E) tends to a con- 
stant for large E; it is not normalizable. This means that, 
for X < Xg, the system does not have a steady state; it 
is ("weakly") non-ergodic and "ages" by evolving into 
deeper and deeper traps |12 1^. The model (|l|) therefore 
has a glass transition at x = Xg. 

With Bouchaud's model, we have a good candidate 
for describing in a relatively simple way the glassy fea- 
tures of SGMs. Its disadvantages for our purposes are: 
(i) The assumption of an exponentially decaying p{E) is 
rather arbitrary in our context. It can be justified in sys- 
tems with "quenched" {i.e., fixed) disorder, such as spin 
glasses, using extreme value statistics (see e.g. 35|), but 



it is not obvious how to extend this argument to SGMs. 
(ii) The exponential form of the activation factor in (|l|) 
was chosen by analogy with thermal activation. But for 
us, X describes effective noise arising from interactions, 
so this analogy is by no means automatic, and functional 
forms other than exponential could also be plausible. In 
essence, we view (i) together with (ii) as a phenomenolog- 
ical way of describing a system with a glass transition. 



B. Incorporating deformation and flow 

To describe deformation and flow, the SGR model |1^] 
incorporates strain degrees of freedom into Bouchaud's 
glass model. A generic SGM is conceptually subdivided 
into a large number of mesoscopic regions, and these 
form the "elements" of the model. By mesoscopic we 
mean that these regions must be (i) small enough for a 
macroscopic piece of material to contain a large num- 
ber of them, allowing us to describe its behavior as an 
average over elements; and (ii) large enough so that de- 
formations on the scale of an element can be described 
by an elastic strain variable. For a single droplet in a 
foam, for example, this would not be possible because of 
its highly non-afflne deformation; in this case, the ele- 
ment size should therefore be at least a few droplet di- 
ameters. The size of the elements is chosen as the unit 
length to avoid cumbersome factors of element volume 
in the expressions below. We emphasize that the sub- 
division into mesoscopic elements is merely a conceptual 
tool for obtaining a suitably coarse-grained description 
of a SGM. The elements should not be thought of as 
sharply-defined physical entities, but rather as somewhat 
diffuse "blobs" of material. Their size simply represents 
a coarse-graining length scale whose order of magnitude 
is fixed by the two requirements (i) and (ii) above. 

We denote by I the local shear strain of an element 



(more generally, the deformation would have to be de- 
scribed by a tensor, but we choose a simple scalar de- 
scription). To see how I evolves as the system is sheared, 
consider first the behavior of a foam or dense emulsion. 
The droplets in an element will initially deform elasti- 
cally from the local equilibrium configuration, giving rise 
to a stored elastic energy (due to surface tension, in this 
example |19|] ). This continues up to a yield point, charac- 
terized by a strain ly, whereupon the droplets rearrange 
to new positions in which they are less deformed, thus 
relaxing stress. The mesoscopic strain I measured from 
the nearest equilibrium position {i.e., the one the element 
would relax to if there were no external stresses) is then 
again zero. As the macroscopic strain 7 is increased, 
I therefore executes a "saw-tooth" kind of motion [p6| . 
Neglecting nonlinearities before yielding, the local shear 
stress is given by kl, with k an elastic constant; the yield 
point defines a maximal elastic energy E = 2^'y • The 
effects of structural disorder are modeled by assuming a 
distribution of such yield energies E, rather than a sin- 
gle value common to all elements. A similar description 
obviously extends to many others of the soft materials 
mentioned above. 

To make the connection to Bouchaud's glass model, 
yield events can be viewed as "hops" out of a trap (or 
potential well), and the yield energy E is thereby iden- 
tified with the trap depth. As before, we assume that 
yields (hops) are activated by interactions between dif- 
ferent elements, resulting in an effective temperature x. 
The activation barrier is now E — ^kP, the difference 
between the typical yield energy and the elastic energy 
already stored in the element. 

For the behavior of elements in between rearrange- 
ments, the simplest assumption is that their strain 
changes along with the macroscopically imposed strain 
7. This means that, yield events apart, the shear rate 
is homogeneous throughout the material; spatial fiuctu- 
ations of the shear rate are neglected in what can be 
viewed as a further mean-field approximation. The SGR 
model therefore applies only to materials which can sup- 
port macroscopically homogeneous flows (at least in the 
range of shear rates of practical interest). In fact, we 
regard this requirement as a working definition of what 
is meant by a "soft" glassy material. A "hard" glassy 
material, on the other hand, might fail by fracture and 
strong strain localization rather than by homogeneous 
flow. Whether a link exists between this distinction and 
the classification of structural glasses into fragile versus 
strong |33| is not clear to us at present. 

While the SGR model assumes a spatially homoge- 
neous strain rate, it does admit inhomogeneities in the 
local strain I and stress a = kl |37|| . These arise because 
different elements generally yield at different times. To 
describe the state of the system at a given time, we there- 
fore now need to know the joint probability of finding 
an element with a yield energy E and a local strain I. 
Within the SGR model p^, this probability evolves in 
time according to 
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To e-(^-5'^-'')/^ P + Tit) p{E)5{l) 



(3) 



The first term on the rhs describes the motion of the 
elements between rearrangements, with a local strain rate 
equal to the macroscopic one, I = 7. The interaction- 
activated yielding of elements (which is assumed to be 
an instantaneous process on the timescales of interest 
to us) is reflected in the second term. The last term 
incorporates two assumptions about the properties of an 
element just after yielding: It is unstrained {I = 0) and 
has a new yield energy E randomly chosen from p{E), 
i.e., uncorrelated with its previous one. Finally, the total 
yielding rate is given by the analog of (||), 

r(() = r„(e-<^-i"">'-) 

= r„/.^.,P,^.,0.--'«->'- (4) 

Eq. (S) tells us how the state of the system, described 
by P{E,l,t), evolves for a given imposed macroscopic 
strain 7(i). What we mainly care about is of course the 
rheological response, i.e., the macroscopic stress. This is 
given by the average of the local stresses 



a{t) = k{l)p = k jdE dl P{E, I, t) I 



(5) 



Eqs. (|[]^) define the SGR model, a minimal model for 
the rheology of SGMs: It incorporates both the "glassy" 
features arising from structural disorder (captured in the 
distribution of yield energies E and local strains I) and 
the "softness" : for large macroscopic strains, the material 
flows because eventually all elements yield. An intuitive 
picture of the dynamics of the SGR model can be ob- 
tained by viewing each element as a "particle" moving 
in a one-dimensional piecewise quadratic potential, with 
noise- induced hops which become increasingly likely near 
the edge of a potential well (see Fig. S). This also shows 
the hysteresis effects associated with yielding: Once a 
hop to a new well has taken place, a finite strain reversal 
is in general needed before a particle will hop back to its 
old well [|§. 

Before moving on to the exact solution of the SGR 
model, we briefly mention some of its limitations. Among 
the most serious of these is the assumption that the noise 
temperature x and the attempt frequency Fq are con- 
stant parameters of the model. In general, they may be 
expected to depend on the imposed shear rate 7, for ex- 
ample, or in fact have their own intrinsic time evolution. 
In particular, it must be born in mind when interpreting 
our results below that the effective noise temperature x is 
not a parameter that we can easily tune from the outside; 
rather, we expect it to be determined self-consistently by 
the interactions in the system. We discuss these points 
in some detail in Sec. VI, where we also speculate on 
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FIG. 2. Potential well picture of the dynamics of the SGR 
model. Note that the relative horizontal displacement of the 
quadratic potential wells is arbitrary; each has its own inde- 
pendent zero for the scale of the local strain I. The solid verti- 
cal bars indicate the energy dissipated in the "hops" (yielding 
events) from 1 to 2 and 3 to 4, respectively. 



the physical origin of the model parameters x and Fq. 
Within the SGR model, the "prior" density of yield en- 
ergies, p{E), is likewise taken to be a constant. This 
implies the assumption that the structure of the mate- 
rial considered is not drastically altered by an imposed 
flow, and excludes effects such as shear-induced crystal- 
lization. 

The SGR model is also essentially a low-frequency 
model. This is due to our assumption that each element 
behaves purely elastically until it yields and a rearrange- 
ment takes place. In reality, the rheological response of 
an element will be more complex. After the application 
of a strain, for example, there may be a fast relaxation 
of the local stress from its instantaneous value, due to 
local relaxation processes. In a foam, for example, these 
might correspond to small shifts of the bubble positions; 
in the language of mode-co upling theory, they could be 
described as /?- relaxations |32|j39| ]. Such local stress re- 
laxation processes are expected to take place much faster 
than actual yielding events, which involve a more drastic 
reorganization of the structure of the material. For fre- 
quencies smaller than the attempt frequency for yielding, 
tj < Fq, they can therefore be neglected. This then im- 
plies that the elastic properties that we ascribe to local 
elements are those that apply once all fast local stress re- 
laxation processes are complete. We have also neglected 
viscous contributions to the local stress; in foams, for ex- 
ample, these are due to the flow of water and surfactant 
caused by the deformation of the elements. In the low 
frequency regime of interest to us, such viscous effects 
are again insignificant (see e.g. uW), whereas at high 
frequencies the model (^-^ would have to be modified 
appropriately to yield sensible predictions. 

Another restriction of the model is the assumption that 
the elastic constant k is the same for all elements. This 
may not be appropriate, for example, for strongly poly- 
disperse materials; we plan to investigate the effects of 
variable k in future work. We have also made the sim- 



plifying assumption that an element is always unstrained 
directly after yielding. Interaction between neighboring 
elements may however frustrate the relaxat ion to the new 
equilibrium state; we discuss briefly in Sec. IV C how this 
feature can be incorporated into the model. 

Finally, the treatment of energy dissipation during 
yield events within the SGR model may also have to be 
refined. This can be seen by expressing the work done 
on the system in the following way: We multiply the 
equation of motion (g) by the elastic energy ^kP of an 
element and integrate over I and E. Integration by parts 
of the 7 term then just gives the stress (ph, hence 



piE) = exp(-£;) 



(7) 
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(6) 



where the averages are over P{E,l,t). The Ihs is the 
rate of energy input into the system. The first term on 
the rhs, which is a complete time differential, describes 
the part of this energy that is stored as elastic energy 
of the elements. The second term, which is always non- 
negative, is the dissipative part. It is just the average 
over all elements of their yielding rate times the energy 
dissipated in a rearrangement, which we read off as ^kP. 
This means that within the model, every rearrangement 
dissipates exactly the elastic energy stored within the 
element when it yields (see Fig. g). 

In general, this is not implausible. But it implies 
that some rearrangements — those of unstrained (I = 0) 
elements — have no dissipation associated with them [Q . 
In reality, however, the local reorganization of a mate- 
rial during any yield event would always be expected to 
dissipate some energy. How much might depend, for ex- 
ample, on the height of the activation barrier for yield- 
ing, E — ifcZ^. The model in its present form does not 
capture such effects; in fact, the yield energies E do not 
feature in the energy balance (o) except through their 
effect on the yielding rates. This exposes a related lim- 
itation of the model: On physical grounds, one would 
expect that elements with a larger yield energy E may 
have a more stable configuration with lower total energy 
(for example, an arrangement of droplets in an emulsion 
with a lower total surface energy). The average value 
of E (which increases during aging, for example p^ , p^ ), 
should then also occur in the energy balance (||) . This is 
not accounted for in the model in its present form. 



III. CONSTITUTIVE EQUATION 

To simplify the following analysis of the model, we 
choose appropriate units for energy and time; a con- 
venient choice is such that Xg = Fq = 1. From the 
definition of the glass transition temperature, this im- 
plies that the density of yield energies has the form 
p{E) = exp[-£:(l + f{E))] with f{E) -^ iov E ^ oo. 
For our numerical investigations below we use the sim- 
plest p{E) of this form, which is purely exponential 



Analytical results, on the other hand, hold for general 
p{E) unless otherwise stated. We eliminate a final pa- 
rameter from the model by setting k = 1; this can always 
be achieved by a rescaling of the stress a and the strain 
variables 7 and I. With this choice of units, it becomes 
clear that the SGR model is in fact rather parsimonious: 
apart from scale factors, its predictions are determined 
by a single parameter, the effective noise temperature 

Note that in our chosen units, typical yield strains 
y^2E/k are of order one. Experimentally, SGMs gen- 
erally have yield stresses of at most a few percent (see 
£.17. [1^,^,Q); the necessary rescaling of our results for 
strain variables should be born in mind when comparing 
to experimental data. For example, a strain rate 7 = 1 in 
our units corresponds to 7 = ^yTo in dimensional units, 
with ly — (xg/kY'^ a typical ("a priori", i.e., sampled 
from p{E)) yield strain. For a specific material, the three 
scale parameters x^, k and Fq of the SGR model could be 
estimated from measurements of a yield strain, a shear 
modulus and a viscosity, for example. 

The derivation of the exact constitutive equation (CE) 
for the SGR model is given App. ^. For simplicity, we 
impose the mild restriction that the initial state is com- 
pletely unstrained, i.e., 7(i = 0) = and 



P{E,l,t^Q) = Po{E)5{l) 



(8) 



Our central result then relates the stress at time t to the 
strain history 7(i') (0 < t' < t) by the CE 

a(i)=7(0Go(^(t,0)) 

+ f dt'T{t')[^{t)~^{t')]G,{Z{t,t')) (9) 



with the yielding rate T(t) determined from 



l = Go(^(t,0))+ / dt' T{t')Gp{Z{t,t')) (10) 

Jo 



Here the functions 



Go{z) = JdEPo{E)eicp(^ 
Gp{z) = I dEp{E)eiip (-ze 



-ze-^/- 



-E/x 



(11) 



describe the purely noise induced decay of the stress. 
This decay is however governed not simply by the time 
interval between a change in macroscopic strain at t' and 
a stress measurement at t, but by an "effective time in- 
terval" z = Z{t, t') given by 

Z(i, t') = 1^* dt" exp { [7(i") - 7(i')]' I2x] (12) 

One reads off that Z{t,t') > t — t'\ the effective time 
interval is always greater than the actual time interval, 



and the more so the larger the changes in strain 7(t") 
from its value at the earlier time t' . This implies a faster 
decay of the stress, and so Z{t, t') can be said to describe 
strain- induced yielding (in other words, shear-thinning). 
In fact, a look at (pUlQ) confirms that all nonlinear ef- 
fects within the model arise from this dependence of the 
effective time interval Z{t,t') on the macroscopic strain 
history 7(t"). 

The CE (|jlO) can be most easily understood by view- 
ing the yielding of elements as a birth-death process: 
Each time an element yields, it "dies" and is "reborn" 
with / = 0. In between such events, its local strain just 
follows the changes in global strain 7(t). If an element 
was last reborn at time t' , its local strain at time t is 
therefore I = ■y{t) — j{t'). Since we set fc = 1, this is 
also its contribution to the stress. The first term on the 
rhs of (H,^0|) is the contribution of elements which have 
"survived" from time to t; they do so with the "sur- 
vival probability" Go{Z{t,Q)). The second term collects 
the contribution from all elements which have yielded at 
least once between time and t, and were last reborn at 
t'. The number of such elements is proportional to the 
rate of "rebirths" at t', i.e., the yielding rate r(t'), and 
the corresponding survival probability Gp{Z{t, t')). Note 
that there are two different survival probabilities here, 
given by Go and Gp, respectively. The difference arises 
from the fact that these probabilities are in fact aver- 
ages over the distribution of yield energies, as expressed 
by (KW- For elements that have survived from t' = 0, 
this distribution is Pq{E), while for elements that have 
yielded at least once, it is p{E). 

The glass y features of the SGR model as discussed 
in Sec. [I A are reflected in the CE (p|pX|), in particu- 
lar in the asymptotic behavior of Gp{z). For the sim- 
ple exponential form (^) of p{E), one easily finds that 



Gp{z) 



— ™i 



asymptotically. As shown in Appendix f 



the same behavior holds for general p{E), in the sense 
that 



hm Gp{z)z^ 

z — >oo 

lim Gp{z)z^ 







(13) 



for any arbitrarily small e > 0. Wc shall refer to this 
property by saying that Gp{z) decays asymptotically as 
z~^ up to "sub-power law factors". Unless otherwise 
specified, all power laws referred to in the following hold 
for general p{E), up to such sub-power law factors. 

Consider now the case where strain-induced yielding 
can be neglected, such that Z{t, t') = t—t'. This is always 
true for sufficiently small strain amplitudes. Below the 
glass transition {x < 1), the time integral /„ dt' Gp{t~t') 
of the response function Gp{Z{t,t')) — Gp(t — t') in (g) 
then diverges in the limit t -^ oo. Compatible with the 
intuitive notion of a glass phase, this means that the sys- 
tem has a very long memory (of the kind that has been 
described as "weak long term memory" Q,Q) and is 
(weakly [n2|) non-ergodic. This can lead to rather intri- 
cate aging behavior, which we plan to explore in future 



work. For the purpose of the prese nt pap er — with the 

we focus on 
These include 
> 1, and the 



exception of a brief discussion in Sec. IV B 



situations where the system is ergodic. 
the regime above the glass transition, x 
case of steady shear flow for all noise temperatures x 
(strain-induced yielding here restores ergodicity even for 
X < 1). In the former case, a choice needs to be made for 
the initial distribution of yield energies. We consider the 
simplest case where this is the equilibrium distribution 
at the given x, 

PoiE) = Peq(S) = reqexp(^/a:)p(^) (14) 

Correspondingly, we write Go{z) = Gcq{z). The function 
Gp{z) is then related to the derivative of G'oq(z) by 

Gp{z) = -T-^'G',^{z) (15) 

with a proportionality constant given by the equilibrium 
yielding rate 

r-qi = fdE p{E) exp{E/x) ^ I dz Gp{z) (16) 



IV. LINEAR RESPONSE 
A. Above the glass transition 

The simplest characterization of the rheological behav- 
ior of the SGR model is through its linear rheology This 
describes the stress response to small shear strain per- 
turbations around the equilibrium state. As such, it is 
well defined (i.e., time-independent) a priori only above 
the glass transition, a; > 1 (see however Sec. IV B). 

To linear order in the applied strain 7(t), the effec- 
tive time interval Z{t,t') — t — t' . In the linear regime, 
all yield events are therefore purely noise-induced rather 
than strain-induced. Correspondingly, the yielding rate 
as determined from (|lO|) is simply r{t) ~ Foq, as can be 



confirmed from eqs. pBUlq ). The expression (y) for the 
stress can then be simplified to the familiar form 



a(i)= [ dt'j{t')Gc 
Jo 



,{t-t') 



(17) 



As expected for an equilibrium situation, the response 
is time-translation invariant p6[ , with Gcq{t) being the 
linear stress response to a unit step strain at i = 0. The 
dynamic modulus is obtained by Fourier transform. 



G*{uj) =iuj I dt e-''^* G^q{t) 
Jo 



lOJT 
iuJT + 1 



(18) 



oq 



This an average over Maxwell modes with relaxation 
times T. For an element with yield energy E, t = 
exp{E/x) is just its average lifetime, i.e., the average 
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FIG. 3. Linear moduli G' (solid line) and G" (dashed) 
vs frequency ui at various noise temperatures x. We only 
show the behavior in the low frequency regime oj < 1, where 
the predictions of the SGR model are expected to be phys- 
ically relevant. The high frequency behavior (predicted as 
G" ~ const, G" ~ i^~^) is not realistic because the model ne- 
glects local viscous effects (among others) which can become 
important in this regime. 



time between rearrangements. The relaxation time spec- 
trum therefore follows from the equilibrium distribu- 
tion of energies, Pcq[E) ex eyip{E / x) p{E) . Because of 
the exponential tail of p{E), it has a power-law tail 
Pcq{T) ^ T~^ (for r 3> 1, up to sub-power law factors). 
As X decreases towards the glass transition, this long- 
time part of the spectrum becomes increasingly domi- 
nant and causes anomalous low frequency behavior of 
the moduli, as shown in Fig. 0: 



G" - uj^ for 3<x, 
G" r^uo for 2 < x. 



■ uj'-^ ^ for 1 < X < 3 

' uj'^-^ for 1 < X < 2 (19) 



For a; > 3 the system is Maxwell-like at low frequencies, 
whereas for 2 < a; < 3 there is an anomalous power law 
in the elastic modulus. Most interesting is the regime 
1 < a; < 2, where G" and G" have constant ratio; both 
vary as uj^^^. Behavior like this is observed in a num- 
ber of soft materials |4|-|7|Jll](] . Moreover, the frequency 
exponent approaches zero as x -^ 1, resulting in es- 
sentially constant values of G" and G", as reported in 
dense emulsions, foams, and onion phases P-p|. Note, 
however, that the ratio G" jG' ~ x — 1 becomes small 
as the glass transition is approached. This increasing 
dominance of the elastic response G" prefigures the on- 
set of a yield stress for x < 1 (discussed below) . It does 
not mean, however, that the loss modulus G" for fixed 
(small) u) always decreases with x\ in fact, it first in- 
creases strongly as x is lowered and only starts decreas- 
ing close to the glass transition (when a; — l~|lnti;|~^). 
The reason for this crossover is that the relaxation time 
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FIG. 4. Linear moduli G' (solid line) and G" (dashed) vs 
frequency a; at a; = 0.9 with energy cutoff -Emax = 10 (thick 
lines) and -Emax = 15 (thin lines) . The loss modulus increases 



as G" 



as the frequency decreases; at very low fre- 



quencies, there is a cross-over to Maxwellian behavior. 



t{{E)\ = exp{{E) /x) corresponding to the mean 



/cq 



equilibrium energy {E) ~ (a; — 1) ^ eventually becomes 
greater than u^^. 



B. Glass phase 

The above linear results only apply above the glass 
transition (a: > 1), where there is a well defined equi- 
librium state around which small perturbations can be 
made. However, if a cutoff -Emax on the yield energies is 
introduced (which is physically reasonable because yield 
strains cannot be arbitrarily large), an equilibrium state 
also exists for x < 1, i.e., below the glass transition. 
(Strictly speaking, with the cutoff imposed there is no 
longer a true glass phase; but if the energy cutoff is large 
enough, its qualitative features are expected to be still 
present.) One then finds for the low frequency behavior 
of the linear moduli: 



G' « const. 



G" 



(20) 



This applies as long as uj is still large compared to the 
cutoff frequency, Wmin = exp(— i^max/s:)- In this fre- 
quency regime, G" therefore increases as uj decreases, 
again in qualitative agreement with some recent exper- 
imental observations ITHIOJ. An example is shown in 
Fig.|. 

The above results relate to the "equilibrium" (pseudo) 
glass phase. The time to reach this equilibrium state is 
expected to be of the order of the inverse of the small- 
est relaxation rate, w~in = exp(i5max/a;)- For large 
Emax, this may be much larger than experimental time 




the glass transition is particularly noteworthy: Whereas 
G"{ui, t) does tend to zero for i ^ cx), it does so extremely 
slowly (as 1/lnt), while at the same time exhibiting an 
almost perfectly "flat" (G" ~ cj^ for small w) frequency 
dependence. Where such an cj-dependence is observed 
experimentally it may well, therefore, correspond to a 
rheological measurement in an out-of-equilibrium, aging 
regime. In order to test this scenario directly, experi- 
ments designed to measure a possible age dependence of 
the linear moduli would be extremely interesting. Such 
experiments would obviously have to be performed on 
systems where other sources of aging (such as coalescence 
in emulsions and foams, evaporation of solvent etc) can 
be excluded; suspensions of microgel beads, hard sphere 
colloids or colloid-polymer mixtures might therefore be 
good candidates. 



FIG. 5. Age-dependence of the dynamic moduli. Shown 
are G' (solid line) and G" (dashed) vs frequency u) a.t x = 1; 
lines of increasing thickness correspond to increasing age of 
the system: t = 10*, 10^, lO'', 10^. Frequencies are restricted 
to the range ujt > 2-n ■ 10, corresponding to a measurement 
of G*{uj,t) over at least ten oscillation periods. Note the 
difference in horizontal and vertical scales; both G' and G" 
have a very "flat" tj-dependence. 



scales, and the non-equilibrium behavior will then be- 
come relevant instead. We give only a brief discussion 
here and refer to a future publication [ [f7[ for more de- 
tails. From the CE (p|,p^), it can be deduced quite gener- 
ally that the stress response to a small oscillatory strain 
7(t) = 7Kexp(iwi) switched on at i = is 

a{t) ^-i^[G*{uj,t)e"^'] 

with a time-dependent dynamic modulus 



C. Frustration 



G*(uj,t)^l 



dt'e~"^'-*-*">T{t')Gp{t-t') (21) 



This modulus is physically measurable only for cot sig- 
nificantly greater than unity, of course, corresponding to 
a measurement over at least a few periods. Here we 
consider the case of an initial distribution of yield en- 
ergies Po{E) = p{E) (hence Gq = Gp), corresponding to 
a "quench" at t = from x ^ cxo to a finite value of 
X. We solve eq. ( [lO|) for the yielding rate T{t) numeri- 
cally and then evaluate G*{uj^t) using (Ell). FigH shows 
the results for a quench to the glass transition {x = 1). 
Not unexpectedly, the frequency dependence of the mod- 
uli follows the same power laws as in the "equilibrium" 
glass discussed above; the amplitude of these, however, 
depends on the "age" t of the system. For x < 1, one 
finds 1 — G*(w,i) ^ {ujtY^'^ Q; this time dependence is 
the same as for the yielding rate T{t) ||l^, and is closely 
related to the aging of the susceptibility in Bouchaud's 
glass model |12]. The behavior of the loss modulus at 



As pointed out in Sec. [IB, the SGR model in its ba- 
sic form (@) assumes that after yielding, each element of 
a SGM relaxes to a completely unstrained state, corre- 
sponding to a local strain of / = 0. This is almost cer- 
tainly an oversimplification: Frustration arising from in- 
teraction of an element with its neighbors will in general 
prevent it from relaxing completely to its new equilib- 
rium state. This leads to a nonzero local strain I directly 
after yielding. This effect can be built into the model 
by replacing the factor 5{l) in (O) by a probability dis- 
tribution q{l] E) of the local strain I after yielding; this 
distribution will in general also depend on the new yield 
energy E of the element. We consider here the case of 
"uniform frustration" , where the strain I after yielding 
has equal probability of taking on any value between —ly 
and Zy, with ly = {2Ey/'^ being the typical yield strain 
associated with the new yield energy. Because values of 
I outside this interval would not make much sense (the 
element would yield again almost immediately) , this sce- 
nario can be regarded as maximally frustrated. 

An exact CE for such a frustrated scenario can still be 
derived, but it is rather more cumbersome than (|9|Jl(]|) 
due to extra integrations over the strain variable I. The 
dynamic moduli, however, can still be worked out fairly 
easily by considering a small perturbation around the 
steady state of (|^) [with S{1) replaced by q{l;E)]. One 
finds 



G*{uj) 



lOJT 



P 



lUJT 
ILOT +1 X {iuJT + 1)2 



where the relaxation times r — exp[(_B — ■^P)/x\ are now 
dependent on both E and Z, and the equilibrium distri- 
bution over which the average is taken is Pcq{E,l) oc 
exp[(_E — ^l'^)/x]p{E)q{l] E). For the uniform frustration 
case, where qil; E) = Q{E - \P)/{d,EY/'^, the dynamic 
moduli are compared with the unfrustrated case in Fig. ra. 
The main effect of frustration is to add a contribution 
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FIG. 6. Effect of frustration. Shown are G' (solid line) and 
G" (dashed) vs frequency uj a.t x = 1.5; results for uniform 
frustration (in bold) are compared with the unfrustrated case 
(thin lines). 



to the relaxation time spectrum near r « 1; this arises 
from elements which have a strain I sa ±/y after yielding 
and therefore yield again with a relaxation rate of order 
unity. Otherwise, however, the main qualitative features 
of the unfrustrated model are preserved; in particular, it 
can be shown that the low frequency power law behav- 
ior (noh remains unchanged. We expect that the same 
will be true for other rheological properties and therefore 
neglect frustration effects in the following. 



V. NONLINEAR RHEOLOGY 

Arguably, the linear rheological behavior described in 
the previous section follows inevitably from the existence 
of a power law distribution of relaxation times. If we 
were only interested in the linear regime, it would be 
simpler just to postulate such a power law. The main 
attraction of the SGR model is, however, that it also 
allows nonlinear rheological effects to be studied in detail. 
It is to these that we now turn. 



A. Steady shear flow^ 

1. Flow curves 

Steady shear flow (7 = const.) is one of the sim- 
plest probes of nonlinear rheological effects. For the SGR 
model, the flow curve (shear stress as a function of shear 
rate) can be calculated either from the long-time limit of 
the CE (|9|,p^), or directly from the steady state solution 
of the equation of motion (0). Either way, one obtains 
for the shear stress 




FIG. 7. Shear stress a vs shear rate 7, for x — 0.25, 0.5, 
. . ., 2.5 (top to bottom on left); x = 1 and 2 are shown in 
bold [Esl. The inset shows the behavior on a linear scale, 
with yield stresses for x < 1 indicated by arrows. 



ct(7) 



l°°dllG,{Z{l)) 
S^dlG,{Z{l)) 



where 



Z{1) = - [ d/'e'"/2- 



1 Jo 



(22) 



(23) 



Eq. (|2^) is just the local strain averaged over its steady 
state distribution, which is proportional to Gp{Z{l)) (for 
I > Q). The resulting stress can easily be evaluated nu- 
merically to give the results in Fig. M. For large shear 
rates 7^1, the shear stress a increases very slowly for 
all X {(J ^ (a;ln7)^'^), corresponding to strong shear 
thinning. More interesting (and more physically rele- 
vant p9[) is the small 7 behavior, where we find three 
regimes: 

(i) For X > 2, the system is Newtonian, a — rj-j, for 
7 — > 0. The viscosity can be derived by noting that in 
this regime, the size of the local strains I that contribute 
significantly to a is proportional to 7. For 7 -^ 0, it 
decreases to zero, and we can approximate Z{1) — l/j, 
giving 

_a _ J^dttGpjt) 



T,^fdEp{E)e 



2E/x 



,E/x 



(r) 



oq 



The viscosity is therefore simply the average of the relax- 
ation time r = eyip{E/x) over the equilibrium distribu- 
tion of energies, Peq{E) = Teqexp{E/x)p{E). From the 
form 77 oc {ex-p{2E/x)) one sees that it diverges at a; = 2, 
i.e., at twice the glass transition temperature. The exis- 
tence of several characteristic temperatures in the SGR 



model is not surprising; in fact, Bouchaud's original glass 
model already has this property [[l3| (which has also been 
discussed in more general contexts, see e.g. |pO[). 

(ii) The divergence of the viscosity for x ^ 2 signals 
the onset of a new flow regime: for 1 < x < 2 one finds 
power law fluid rather than Newtonian behavior. The 
power law exponent can be derived as follows: The steady 
shear stress (P2) is the ratio of the integrals 



In[i) 



dl rGp{Z{l)) 



for n — 1 and n = 0. By techniques very similar to those 
used in App. g, one derives that in the small 7 limit, 
/„ scales as 7"+^ for x > n + 1; for lower x, there is 
an additional contribution scaling as 7^ up to sub-power 
law factors (see App. O). The dominant contribution to 
a for small 7 in the regime 1 < a; < 2 therefore scales as 
0- r^ 7^~^, again up to sub-power law factors. The power 
law fluid exponent therefore decreases linearly, from a 
value of one for a; — 2 to zero at the glass transition 
x = 1. 

(iii) For x < 1, the system shows a yield stress: (7(7 -^ 
0) = (Ty > 0. This can again be understood from the 
scaling of /i and Iq: the dominant small 7 contributions 
to both scale as 7^ for x < 1, giving a finite ratio fXy = 
h/Io in the limit 7 — + 0. For general p{E) there are 
subtleties due to sub-power law corrections here, which 
are discussed in App. O. Here we focus on the simplest 
case (0) of exponential p{E), where such corrections are 
absent. Using the scaling of /i and /q, we can then write 
the shear stress for small 7 as 




FIG. 8. Yield stress CTv ss a function of x. 
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(24) 



Beyond yield, the stress therefore again increases as a 



FIG. 9. Effective power law exponent dln(cr — cry)/dln7 vs 



power law of the shear rate, cr — CTy ex 7 . For ex- 7 in the glass phase (left, yield stress o-y > 0, a: = 0.1, 0.2 



ponential p{E), the yield stress itself can be calculated 
explicitly: In order to have dy > 0, the values of I that 



contribute to the shear stress (22) must remain finite for 
7^0. But then for any fixed I, Z{1) -^ 00. We can 

x\z^ 



therefore use the asymptotic form Gp{z) 



— 7*1 y ^ 



in (22) 



givmg 



Cdll[Z{l)]- 
Io"dl[Z{l)]- 



(25) 



The factors 7^ (from the definition (|2^) of Z{1)) in nu- 
merator and denominator have canceled jinaking the re- 
sult independent of 7 as required. Fig. |8| shows the re- 
sulting yield stress as a function of a:; it has a linear onset 
near the glass transition, cry '^ 1 — x. 

To summarize, the behavior of the SGR model in 
regimes (ii) and (iii) matches respectively the power-law 
fluid pl-pi and Herschel-Bulkeley 0,0| scenarios as used 
to fit the nonlinear rheology of pastes, emulsions, slurries, 
etc. In regime (ii), the power law exponent is simply x—1, 
X being the effective (noise) temperature; in regime (iii) 



. . . 0.9 from top to bottom) and in the power law fluid regime 
(right, ay — 0, X = 1.1, 1.2 . . . 1.9 from bottom to top). 



and for exponential p{E), it is 1 — a; (see App. ^ for a 
discussion of the general case). Numerical data for the 
effective exponent dln{a — ay)/dlnj in Fig. |9| are com- 
patible with this, although the exponent only approaches 
its limiting value very slowly as 7 — > for x near the 
boundaries of the power law regime, x — 1 and 2. 

A natural question to ask is of course how the existence 
of a yield stress in the glass phase affects the linear mod- 
uli, i.e., the response to small strains. This is a highly 
nontrivial issue due to the non-ergodicity of the glass 
phase and the corresponding aging behavior. In particu- 
lar, the answer will depend to a significant degree on the 
strain history of the material. We therefore leave this 
point for future, more detailed study |47[. 
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2. Flow interrupts aging 



We saw above that there is a steady state regime for 
any value of x in the presence of steady shear flow. On 
the other hand, the discussion in Sees. HI and IV B 
showed that in the absence of flow, the system has no 
steady state in the glass phase (a; < 1) and instead ex- 
hibits aging behavior. The difference between the two 
cases can be seen more clearly by considering the distri- 
bution of yield energies, P{E). Without flow, one obtains 
a Boltzmann distribution P{E) ex p{E)eKp{E/x) up to 
(for x < 1) a "soft" cutoff which shifts to higher and 
higher energies as the system ages |1^. This cutoff, and 
hence the most long-lived traps visited (which have a life- 
time comparable to the age of the system), dominate the 
aging behavior 11^]. In the presence of flow, on the other 
hand, there is a finite steady state value for this cutoff; 
one finds 

P{E) ex p{E) e^/^ for E -^ xlni^-^x^/"^) 

P{E) ex p{E) E^l'^ for E ~> xX^ici^^x^''^) (26) 

(only the second regime exists for 7 > x^l"^). The ex- 
istence of these two regimes can be explained as fol- 
lows: Assume the yielding of an element is noise-induced. 
Its typical lifetime is then exp(i?/a:), during which it 
is strained by 7exp(_E/a;). The assumption of noise- 
induced yielding is self-consistent if this amount of strain 
does not significantly enhance the probability of yielding, 
i.e., if [7exp(i?/x)]^/a; <C 1. This is the low E regime 
in (Eq), which gives a Boltzmann form for the yield en- 
ergy distribution as expected for noise-induced yielding. 
In the opposite regime, yielding is primarily strain in- 
duced, and the time for an element to yield is of the 
order of Zy/7 = (2_B)^/^/7 (rather than exp(_E/x)). In- 
tuitively, we see that flow prevents elements from getting 
stuck in progressively deeper traps and so truncates the 
aging process after a finite time. We can therefore say 
that "flow interrupts aging" pj]. 



3. Cox-Merz rule 

A popular way of rationalizing flow curves is by relat- 
ing them to the linear rheology via the heuristic Cox- 
Merz rule [ pl| . This rule equates the "dynamic viscos- 
ity" 77* (w) = \G*{iS)\/ui with the steady shear viscosity 
77(7) = <^{i)li when evaluated at 7 = lj. The ratio 
bJrj{'^ = uj) / \G* {ijj)\ is therefore equal to unity if the Cox- 
Merz rule is obeyed perfectly. Using our previous re- 
sults, we can easily verify whether this is the case in the 
SGR model. From Fig. |l^, we see that in the Newtonian 
regime x > 2, the Cox-Merz rule is obeyed reasonably 
well for frequencies w < 1; for cj — > 0, it holds exactly 
as expected (recall that 77(7) = (t), while from (19), 
G*{uj — » 0) = zcij(r)). In the power-law fluid regime 
1 < X < 2, on the other hand, the Cox-Merz rule is seen 
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FIG. 10. Cox-Merz ratio 0^77(7 = lj)/\G*{uj)\ as a function 
of uj for noise temperatures x = 1, 1.2, . . . , 1.8, 2 (bold), 2.5, 
3 (bottom to top). 



to be less reliable and is not obeyed exactly even in the 
zero frequency limit. At the glass transition (x — > 1), 
it fails rather dramatically: In this limit, \G*{Ld)\ = 1 
and so the Cox-Merz rule predicts a shear rate indepen- 
dent shear stress (7(7) — 777(7) — 1, whereas in fact (7(7) 
decreases to zero for 7^0. 



4.. Dissipation under steady shear 

Finally, in conclusion of this section on steady shear 
flow, we calculate the distribution of energies dissipated 
in yield events. This distribution may provide a use- 
ful link to computer simulations of steady shear flow of 
foams, for example, where it is often easy to monitor dis- 
continuous drops in the total energy of the system and 
determine their distribution |23|| . The correspondence is, 
however, not exact. Our mean-field model treats all yield 
events as uncorrelated with each other, both in time and 
space. In reality, such correlations will of course exist. 
In fact, several events may occur simultaneously, at least 
within the time resolution of a simulation or experiment. 
The observed drop in total energy would then have to 
be decomposed into the contributions from the individ- 
ual events to allow a direct comparison with our model. 
This is only possible if the events are sufficiently localized 
(spatially) to make such a decomposition meaningful. In 
foams and emulsions, there is evidence that this may in- 
deed be the case @|3|j42]J5|-|§. 

We earlier derived the energy balance equation (|6|) and 
deduced from it that, within the model, each yield event 



dissipates the elastic energy /S.E — ^P stored in the el- 
ement just prior to yielding. The probability of observ- 
ing a yield event with energy dissipation AE is therefore 



11 




FIG. 11. Distribution P{AE) of energies AE dissipated in 
yield events under steady flow, for x = 1.5 and 7 — 10^*, 
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given by 



1 (bottom to top at AE = 1) 



P{AE) ^ - 



dEdl P{E, I) e-(^-5'')/^ ^J'a£; _ h^ 



The steady state distribution P{E, I) of yield energies 
and local strains for a given shear rate 7 and noise tem- 
perature X can easily be deduced from (0). After some 
algebra, the result can be put into the simple form 

P{AE) dAE = -4 GJZa)) dl 
ol 

Fig. |ll| shows the resulting P{AE) for exponential p{E). 
Larger shear rates 7 are seen to lead to an increas- 
ing dominance of "large" yield events (which dissipate 
a lot of energy). This is intuitively reasonable: the 
larger 7, the larger the typical strains of elements when 
they yield. The functional dependence of P{AE) on 
AE is surprisingly simple. An initial power law decay 
P{AE) ^ AE^^/"^ crosses over for AE « 7^ into a sec- 
ond power law regime P{AE) ^ AE^^^^I"^ . This is 
cut off exponentially for values of AE around unity [^ . 
The exponential tail for very large dissipated energies is 
P{AE) ~ exp(— A£') independently of x. This asymp- 
totic behavior is the same as for the prior density of yield 
energies, p{E) ~ exp(— £■); measurements of P{AE) for 
large AE could therefore yield valuable information on 
piE). 

These results for P{AE) also help to understand the 
small 7 scaling of the energy dissi pation rate a-y — 
r (AE). From the results of Sec. V A , we know that this 
is 7^ in the Newtonian regime a; > 2, 7^ in the power law 
fluid range 1 < a; < 2, and 7 in the yield stress regime 
X < 1. (The limit 7 ^ is always understood here and 
in the following.) The form of P{AE) suggests to decom- 
pose the dissipation into its contributions from "small" 
{AE = 0(7^)) and "large" (A£' = 0(1)) dissipation 



events. Each of these two classes makes a contribution 
to 177 which is the fraction of elements in the class, times 
the average yielding rate in the class, times the average 
energy dissipated. Hence, in obvious notation, 

ai = PsTsAEs+PiTiAEi 

One then easily confirms the following scalings. The av- 
erage dissipated energies are obviously given by AEs = 
0(7^) and AEi — 0{1). The average yielding rate for the 
small, noise induced events is independent of shear rate, 
Ts = 0(7°); while for the large, shear induced events it 
is F; = 0(7). Finally, for the fractions of small and large 
elements, one finds that above the glass transition, almost 
all elements have small strains I = 0(7), corresponding 
to AE = 0(7^); hence Pg = 0(1). Large strains, on the 
other hand, occur with a probability Pi = 0(7^"^) which 
becomes vanishingly small for small shear rates. Below 
the glass transition, the situation is reversed: Pi = 0(1), 
while Pg = 0(7^^^). Putting everything together, one 
has: 

(i) In the Newtonian regime {x > 2), dissipation is 
dominated by small, noise induced events, and is there- 
fore of 0(7^). 

(ii) In the power law fluid range (1 < a:; < 2), a 
vanishingly small number of elements has large strains, 
but these dominate the dissipation 177 — PiTiAEi — 
0(7'^~^)0(7) — 0['y^). As the glass transition is ap- 
proached, the fraction of large elements and hence the 
dissipation increases. 

(iii) In the yield stress regime, most elements have large 
strains, giving a dissipation rate (77 — 0(7) which simply 
scales with the shear rate. 

With the same approach, one can also analyse the total 
yielding rate F — PsTs + PiTi. Small events always dom- 
inate, and F therefore scales with 7 in the same way as 
Ps ■ This is true even in the non-Newtonian flow regimes 
{x < 2), where the contribution of these elements to the 
total dissipation rate is negligible. 

The distribution of total energy drops AStot due to re- 
arrangements has been monitored in recent simulations 
of steady shear flow of two-dimensional foam, based on 
a "soft-sphere model" [g2|,g3|. It was found to exhibit a 
power law P(AEtot) ^ ^^tot with an exponent 1/ « 0.7, 
with an exponential cutoff for large energy drops. More 
recent simulations using the same model suggest that, 
when AStot is normalized by the average elastic energy 
per foam bubble, the form of P(A£'tot) is largely insen- 
sitive to variations in shear rate 7. Decreasing the gas 
volume fraction (/), on the other hand, moves the (normal- 
ized) cutoff to larger energies, suggesting a possible di- 
vergence near the rigidity loss transition at (j) ~ 0-64 p8[ . 
Simulations using a "vertex model" , on the other hand, 
gave P{AEtot) '^ ^-E-tot with no system-size indepen- 
dent cutoff for large Ai^tot |2l[l . It is unclear how these re- 
sults can be reconciled; neither, however, is fully compat- 
ible with the predictions of the SGR model for P{AE). 
At this point, we do not know whether this disagreement 
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is due to the difference between IS.E (dissipation in a sin- 
gle yield event) and Ai?tot (total dissipation in a number 
of simultaneous yield events), or whether it points to a 
more fundamental shortcoming of the SGR model such 
as neglect of spatial or temporal correlations. 



B. Shear startup 

If a shear flow is started up at t = 0, such that 
7(t) = 7^ for i > 0, then (7(7) as given by the flow curve 
is the asymptotic, steady state value of the stress for 
t — *■ 00. We now consider the transient behavior a{t) for 
finite t. This depends on the initial state of the system 
at i = 0; here we consider only the case where this ini- 
tial state is the equilibrium state (|lj) at the given value 
of X. This restricts our discussion to the regime above 
the glass transition, x > 1, where such an equilibrium 
state exists [^. Solving the CE (P,[lO|) numerically, we 
can find the stress cr as a function of time t or, alterna- 
tively, strain 7. Fig. |l2| shows exemplary results. The 
initial behavior under shear startup is found to be elas- 
tic in all cases, a — ^. (This can in fact be deduced 
directly by expanding (0) to first order in t and noting 
that G[){Z{t,Q)) = 1-1- 0{t) while the contribution from 
the integral is of 0{t^).) Asymptotically, on the other 
hand, the stress approaches the steady-state (flow curve) 
value (7(7). However, the model predicts that it does 
not necessarily do so in a monotonic way. Instead, the 
stress can "overshoot"; within the model, this effect is 
most pronounced near the glass transition {x w 1). Such 
overshoot effects have been observed experimentally in, 
for example, foam flow [p[. The tendency towards large 
overshoots for a; ^ 1 agrees with our results for the lin- 
ear moduli and flow curves: As the glass transition is 
approached, the behavior of the system becomes predom- 
inantly elastic; the stress can therefore increase to larger 
values in shear startup before the material (as a whole) 
yields and starts to flow. 



C. Large step strains 

As a further probe of the nonlinear rheological behav- 
ior predicted by the SGR model, we now consider large 
(single and double) step strains. Again, we do not dis- 
cuss aging effects here and therefore limit ourselves to the 
regime x > 1 with the equilibrium initial condition (14). 

The case of a single step strain {'^{t) = 76(i), with 
0(i) = 1 for t > and zero otherwise) is particularly 



simple. The integral over t' in the CE (O) is then identi- 
cally zero, giving a stress response of 

ait) - 7Go(^(t, 0)) - 7Gcq (e^'/'^i) (27) 

Comparing with the response (O) in the linear regime, 
the effect of nonlinearity is to speed up all relaxation pro- 
cesses by a factor exp(7^/2x). It is easy to see why this 




FIG. 12. Stress a vs strain 7 for shear startup at effective 
temperature x = 1.5. The shear rate 7 = 0.001, 0.002, 0.005, 
0.01, 0.02, 0.05, 0.1 increases from bottom to top. 



is the case. Because we are starting from an unstrained 
equilibrium configuration, each element initially has I = 
and a yielding rate exp(— S/x). Directly after the strain 
is applied, it therefore has local strain I = 7; this in- 
creases its relaxation rate to cxp[— (_E — i^l^)/^]! i-e., by 
the same factor exp(7^/2x) for all elements. Fig. Il^ il- 
lustrates this effect of strain nonlinearity; note that the 
stress for large step strains can decay to small values 
faster than for small strains, due to the strain-induced 
speed-up of all relaxation processes. 

Interestingly, the instantaneous response is always elas- 
tic and not affected by nonlinear effects: a{t = 0+) = 7 
for all 7. It is easily shown from the CE (HJIQ) that 
this is a general feature of the SGR model; whenever 
the macroscopic strain 7(i) changes discontinuously by 
A7, the stress ait) changes by the same amount. We 
also note that the stress response (E^) cannot be fac- 
torized into time and strain dependence. However, for 
the particular case of exponential p{E) and long times 
exp(7^/2x)i ^ 1, such a factorization does exist due to 



the asymptotic behavior of Geo, G'ca(z) 
follows from Gp{z) ^ ; 



(This 



and (p_5|).) One then has 



ait) - 7ft,(7)G'oq(t) h{^) = exp 



-\{l-^-')l' 



The product 7/1(7) tends to zero as 7 increases, corre- 
sponding to a pronounced shear-thinning effect. 

By applying two (large) step strains in sequence, one 
can further probe the nonlinear response of the SGR 
model. Let 71 and 72 be the amplitudes of the two 
strains. If the first strain is applied at t = and the sec- 
ond one at t = At, then 7(t) = 7iB(t)+ j2&{t - At). It 
is straightforward to solve the CE (p|jlC|) numerically for 
t > At. Fig. 113 exemplifies the results for the two cases 
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FIG. 13. Stress response to step strains of amplitude 7 = 1, 
2, 3, at noise temperature x — 1.5. 



where the strains are either equal or of equal magnitude 
but opposite sign. In the first case, and more generally 
when 7i72 > 0, the second step strain speeds up the 
stress relaxation (by a factor exp{[(7i + 72)^ — 7^]/2a;} 
for small At). Therefore, even though the stress is in- 
creased momentarily when the second strain is applied, 
it can actually relax back to zero more quickly than in 
the absence of this strain. In the second case (7172 < 0), 
the second step strain can to some degree reverse the 
speed-up from the first step strain. A particularly sim- 
ple form of the resulting stress response is obtained for 
71 = — 72 = 7 and small At: 



a{t > At) = -7 
X G 



1 - Geq ( e^ /'"At 



(e^'/2-(i-At)) 



This can be understood by noting that the stress for 
t > At is due entirely to elements which have yielded be- 
tween the application of the first and the second strain; 
all other elements have simply followed the two changes 
of macroscopic strain and are therefore back to their un- 
strained state / = after the second strain. The factor 
in squared brackets just gives the fraction of such ele- 
ments. The time dependence of the ensuing stress re- 
laxation is determined by Gp rather than Ggq because 
elements that have yielded were "reborn" with yield en- 
ergies sampled from p{E). These elements — which have 
"forgotten" about the first step strain — also receive a 
speed-up of their relaxation by the second strain. 

The above results can be compared to the predictions 
of the empirical BKZ (Bernstein, Kearseley, Zapas) equa- 
tion Pol . This relation approximates the stress response 
to an arbitrary strain history in terms of the response 
a{t) = (j)(t, 7) to a step strain of size 7 at time t = 0: 
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FIG. 14. Stress response to two large step strains of (a) 
equal (71 = 72 = 2) and (b) opposite (71 = —72 = 2) sign, 
applied at times i = and t = At = 0.1, 0.5, 1, 2, 5, respec- 
tively. Noise temperature x = 1.5. 
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CTBKzit)= [ dt' J^4>{t-t\l) 

For two step strains, this gives, for t > At 



7 = 1 in our units therefore corresponds to a real strain 
of generaUy at most a few percent. 

We consider only the ergodic regime a; > 1; we also 
ignore transient behavior caused by start-up of the oscil- 
latory strain. In the steady state, we can write the stress 



o^EKzli) = (t>it, 7i + 72) - 0(i, 72) + 4>{t - At, 72) (28) response to an oscillatory strain j{t) = j^e^' 



In our case, <p{t,j) is given by (p7|), and the BKZ pre- 
diction is plotted in Fig. |l^ along with the exact results. 
One finds that for the SGR model, the BKZ equation is 
at best approximate, at worst qualitatively wrong. This 
is most easily seen in the size of the stress jump at i = At; 
the BKZ equation predicts 



0(0+, 72) + [0(At,7i +72) - <^(At,7i) 



<^(Ai,72)] 

(29) 



Because (/)(0''",7) = 7 within the SGR model, the term 
in square brackets is the deviation from the true value, 
which is 72. For 71 = —72, the BKZ prediction for the 
stress jump is exact because (/'(t, 7) = — (/'(t, —7); in this 
case (Fig. |ljb), it also works reasonably well for the sub- 
sequent stress relaxation. In the general case, however, 
it is unreliable; Fig. [ija shows that it can in fact even 
predict the wrong sign for the stress jump. 

Finally, we note that a failure of the BKZ equation 
has also been observed in double step strain experiments 
on polymeric liquids | |61| . There, however, the most pro- 
nounced deviations occur for successive step strains of op- 
posite sign rather than, as in the SGR model, for strains 
of the same sign. This can be understood on the basis 
of the different kinds of nonlinearities in the two cases. 
Roughly speaking, in the polymer case the BKZ equation 
fails because it neglects memory of the shape of the tube 
in which a given polymer molecule reptates |£l|,^ . Such 
memory effects are strongest for strain reversal, which 
can bring the tube back to a conformation close to its 
original shape. In the SGR model, on the other hand, 
the BKZ equation fails because it does not adequately 
represent the effects of the strain history on the stress re- 
laxation rates in the material. Such effects are strongest 
when an applied strain compounds an earlier speed-up of 
relaxation processes, i.e., for double step strains of the 
same sign. 



D. Large oscillatory strains 

1. Dynamic moduli 

As a final example of nonlinear rheological behavior, 
we consider the case of large oscillatory strains. We re- 
mind the reader at this point that we have chosen units 
in which typical local yield strains are of order unity (see 
Sec. III). To transform to experimentally relevant quan- 
tities, all strain values have to be multiplied by a typical 
yield strain ly of the SGM under consideration. A strain 



a{t) =7SR[G*(w,7)e" 



Aa{t) 



(30) 



where Aa(t) contains the contributions from all higher 
harmonics. This defines an amplitude dependent dy- 
namic modulus G*{ll>,"/); the relative root-mean-square 
size of the stress contributions from higher harmonics is 
measured by the residual r, defined as 



(31) 



Idt[Aa{t)Y 
Jdta^t) 



The determination of G* and r from the CE (|9UlC 
presents no conceptual difficulties, but is somewhat non- 
trivial numerically (see App. ^ for details). The solution 
yields in fact not just G* and r, but the whole "wave- 
form" of the stress response (T{t). Fig. |5|a shows how 
the response becomes more and more non-sinusoidal as 
the strain amplitude is increased. The stress amplitude 
first increases linearly with 7, then drops slightly as the 
system crosses over from elastic to liquid-like behavior, 
and finally rises again slowly as the typical shear rate 
^Lu of the (now essentially liquefied) material increases. 
Plotting "f{t) and a{t) in a parametric stress-strain plot 
(Fig. [l5|b), one finds a hysteresis loop for large ampli- 
tudes, with stress overshoots near the points where the 
strain rate reverses its sign. 

Consider now the resulting nonlinear modulus G*. 
Fig. [1^ shows an example of a "strain sweep" : The mod- 
uli G'and G" and the residual r are plotted as a function 
of strain amplitude for different frequencies lu. The am- 
plitude dependence of G" is particularly noteworthy: As 
7 increases, G" first increases, but then passes through a 
maximum and subsequently decreases again. This is in 
qualitative agreement with recent measurements of non- 
linear dynamic moduli in, for example, dense emulsions 
and colloidal glasses |0,|lO|,^3|J6J] . The maximum in G" 
is most pronounced near the glass transition x = I; for 
higher noise temperatures, it decreases and disappears 
altogether around x = 2. This is compatible with the 
following coarse estimate of the decay of G" beyond the 
maximum: For sufficiently large strain amplitudes 7, the 
system is expected to flow essentially all the time. If 
the shear rate 7 changes sufficiently slowly {uj <C 1), the 
stress can be approximated as following "adiabatically" 
the instantaneous shear rate: a(t) w cr{'y(t)) with (7(7) 
the steady shear flow curve. For 1 < a; < 2 and suffi- 



ciently small shear rates joj, we know from Sec. V A that 

this relationship is a power law, ct{j) ^ 7^"^. Hence 

a{t) ^ {jajsmajt)^ 

G" ~ 7^-2^ For X 

decay for large 7 (as long as the condition juj <C 1 is 



r 

^ which leads to a 7 dependence of 
-> 2, G" should therefore no longer 
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FIG. 16. Strain sweep; Nonlinear moduli G', G" and resid- 
ual r as a function of strain amplitude 7. Noise temperature 
X = 1.1; lines of increasing thickness correspond to cij = 0.001, 
0.01, 0.1. Recall that 7 is rescaled by a typical local yield 
strain; 7 = 1 therefore corresponds to a real strain of at most 
a few percent. 




FIG. 15. (a) Stress response ait) for oscillatory strain 
7(i) = 'ycos(ut), for frequency u — 0.01 and effective tem- 
perature X — 1.1. Initially, the response is almost perfectly 
elastic; as the strain amplitude increases (curves are shown 
for 7 = 0.1, 0.5, 1, 2, 3, 5), the zero crossings of a{t) move 
to the left, corresponding to progressively liquid-like behav- 
ior (strain lagging behind stress), (b) Parametric plots of 
stress a{t) vs strain 7(t), for same parameter values as in (a); 
7 = 1.5 is also shown. 



obeyed) , in agreement with our observation that its max- 
imum with respect to 7 disappears around this value 
of X. The estimate G" ~ 7^-^ is roughly compatible 
with our numerical data, but a precise verification of this 
power law is difficult (due to severe numerical problems 
for 7 > 20). Note that within the same approximation, 
G' would be estimated to be identically zero, which is of 
course unphysical. Instead, we expect it to decay to zero 
faster than G" as 7 increases, and this is indeed what our 
numerical data show. 



8. Size of linear regime 

The above results allow us to determine the size of the 
linear regime for oscillatory rheological measurements, 
i.e., the largest strain amplitude 7c for which the mea- 
sured values of G" and G" represent the linear response of 
the system. An important first observation that can be 
made on the basis of Fig. 11^ is that the size of the residual 
r is not in general sufficient to determine whether one is 
in the linear regime or not. For example, for strain am- 
plitude 7 = 1.5 at X — 1.1 and lo — 0.1, r is only around 
2.5% even though the value of G" is already twice as 
large as in the linear regime. The a{t) vs "f(t) plot in 
Fig. [l5| b also demonstrates this: for 7 — 1.5, the curve 
still looks almost perfectly elliptical, suggesting linear 
response, while its axis ratio is actually quite different 
from the one in the linear regime. Closer to the glass 
transition, this effect becomes even more pronounced. It 
suggests strongly that whenever the dynamic moduli of 
SGMs are measured, an explicit strain sweep is needed 
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FIG. 17. Frequency dependence of (nonlinear) dynamic 
moduli G'{u),'y) (solid lines) and G"{Lu,"f) (dashed) measured 
at constant finite strain amplitude 7. Noise temperature 
X = 1.001; increasing values of 7 = 0, 1, 2, 3 correspond 
to increasing line thickness. Recall that 7 is rescaled by a 
typical local yield strain; 7 = 1 therefore corresponds to a 
real strain of at most a few percent. The loss modulus G" 
increases strongly with 7, whereas G" varies much less (the 
curves for 7 = and 7 = 1 cannot even be distinguished on 
the scale of the plot). 



to determine whether measurements are actually taken 
in the linear regime. 

If concerns about nonlinear effects are disregarded, an 
experimentally convenient procedure is to measure the 
dynamic moduli at fixed strain amplitude 7 (while vary- 
ing the frequency wV Some numerical results for this case 
are shown in Fig. O. Again, the most interesting behav- 
ior occurs near the glass transition. There, we observe 
that only relatively minor differences in the amplitude of 
the imposed strain can lead to large changes in the mea- 
sured values of G" (whereas G" is affected less strongly). 
This emphasizes again that extreme caution needs to be 
taken in experiments designed to determine the dynamic 
moduli of soft glassy materials; in particular, it needs 
to be born in mind that the loss modulus can easily be 
over-estimated due to undetected nonlinear effects. 

Finally, the actual size of the linear regime itself is also 
of interest. We choose as a working definition of the linear 
regime the strain amplitude 7c at which cither G" or G" 
first deviate by 10% from their values in the limit 7-^0. 
(This implies similar maximum relative deviations for 
\G*\ and the loss tangent tan 5 = G"/G'.) Fig. ^ shows 
7c((^) for several noise temperatures x. Several general 
trends can clearly be read off. First, in the low frequency 
regime, the size of the linear regime decreases as the glass 
transition is approached. This is intuitively reasonable as 
one expects nonlinearities to become stronger near the 
glass transition [p5|. Note, however, that 7c does not de- 



FIG. 18. Size of linear regime 7c vs lo for x = 1.001, 1.5, 2, 
. . . , 5 (bottom to top on left). Close to the glass transition, 
deviations from linearity first show up in G" , which therefore 
determines 7^ (dashed line); for larger x, the linear regime is 
limited by deviations in G' (solid lines). Recall that 7c, like 
all strain variables, is rescaled by a typical local yield strain; 
7c = 1 therefore corresponds to a real strain of at most a few 
percent. 



crease to zero at the glass transition; it tends to a finite 
value of order unity which by our choice of units corre- 
sponds to the typical (a priori) yield stress of local ele- 
ments. The frequency dependence of 7c(cl;) also changes 
as one moves away from the glass transition: Initially 
(for X K, 1), 7c is essentially independent of uj and does 
remain so until around a: = 3 (although its absolute value 
increases); for yet higher noise temperatures, one finds a 
crossover to a 7c ~ i^~^ dependence. The latter corre- 
sponds to the "naive" criterion that the typical shear rate 
^Lo needs to be smaller than typical relaxation rates (of 
order unity away from the glass transition) in order for 
the imposed strain not to create nonlinear effects. The 
predicted w-independence of 7c near the glass transition 
should be easy to verify experimentally. 



VI. INTERPRETATION OF MODEL 
PARAMETERS 

As has been demonstrated above, the SGR model cap- 
tures important rheological features that have been ob- 
served in a large number of experiments, at least in the 
region around the "glass transition" of the model. Us- 
ing a mean-field (one element) picture, it is also simple 
enough to be generic. However, a significant challenge 
that remains is the interpretation of the model parame- 
ters, namely, the "effective noise temperature" x and the 
"attempt frequency" Fq. To tackle these questions, we 
should really start from a more comprehensive model for 
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the coupled nonlinear dynamics of the "elements" of a 
SGM and then derive the SGR model within some ap- 
proximation scheme. At present, we do not know how to 
do this, and the following discussion will therefore have 
to remain rather speculative. 



A. Effective noise temperature x 

We can interpret the activation factor exp[— (£' — 
\kl'^)/x\ in the equation of motion (^ of the SGR model 
as the probability that (within a given time interval of 
order l/Fg) a given element yields due to a "kick" from 
a rearrangement (yield event) elsewhere in the material. 
Therefore x is the typical activation energy available from 
such kicks. But while kicks can cause rearrangements, 
they also arise from rearrangements (whose effects, due 
to interactions, propagate through the material). So 
there is no separate energy scale for kicks: Their energy 
must of the order of the energies released in rearrange- 
ments, i.e., of the order of typical yield energies E. In our 
units, this means that x should be of order unity. Note 
that this is far bigger than what we would estimate if 
X represented true thermal activation. For example, the 
activation barrier for the simplest local rearrangement 
in a foam (a Tl or neighbor-switching process) is of the 
order of the surface energy of a single droplet; this sets 
our basic scale for the yield energies E. Using typical 
values for the surface tension and a droplet radius of the 
order of 1/im or greater, we find E > lO^k^T. In our 
units E = 0(1), so thermal activation would correspond 
to extremely small values oi x = k^T < 10~^. 

We now argue that x may not only be of order one, but 
in fact close to one generically. Consider first a steady 
shear experiment. The rheological properties of a sample 
freshly loaded into a rheometer are usually not repro- 
ducible; they become so only after a period of shearing 
to eliminate memory of the loading procedure. In the 
process of loading one expects a large degree of disorder 
to be introduced, corresponding to a high noise temper- 
ature a; ^ 1. As the sample approaches the steady state, 
the flow will (in many cases) tend to eliminate much of 
this disorder [g6| so that x will decrease. But, as this 
occurs, the noise-activated processes will slow down; as 
x — > 1, they may become negligible. Assuming that, in 
their absence, the disorder cannot be reduced further, x 
is then "pinned" at a steady-state value at or close to 
the glass transition. This scenario, although extremely 
speculative, is strongly reminiscent of the "marginal dy- 
namics" seen in some mean-field spin glass models. In the 
spherical p-spin glass, for example, one finds that after a 
quench from T = cx) to any temperature < T < Tg be- 
low the (dynamical) glass transition temperature Tg, the 
system is dynamically arrested in regions of phase space 
characteristic of Tg itself, rather than the true tempera- 
ture T ||,||. 

There remain several ambiguities within this picture. 



for example whether the steady state value of x should 
depend on 7; if it does so strongly, our results for steady 
flow curves will of course be changed. If a steady flow 
is stopped and a linear viscoelastic measurement per- 
formed, the results should presumably pertain to the x 
characterizing the preceding steady flow (assuming that 
X reflects structure only). But unless the strain ampli- 
tude is extremely small the x-value obtained in steady 
state could be affected by the oscillatory flow itself. This 
might allow "flat" moduh G*[ijj) {x « 1) to be found 
alongside a nonzero yield stress with power law flow ex- 
ponent around 1/2 {x w 1/2) @,||,|3- 

Experimentally, the above ideas concerning the time 
evolution of x in steady flows could be tested in systems 
which can be prepared in both low- and high-disorder 
states, such as onion phases pq |: Strain induced order- 
ing starting from an initial x well below or above Xg = 1 
should drive the system towards x — Q ov x k 1^ respec- 
tively, leading to different rheological characteristics. 

Theoretically, the minimal extension to the SGR model 
that would be needed to substantiate the above scenario 
would be to allow x to evolve in time. We do not know at 
present how to deduce the correct form of this evolution 
in a principled way from some underlying microscopic 
dynamics. However, one possibility is to couple x to the 
number of rearrangements in the material, i.e., the yield- 
ing rate V. Indeed, suppose we view V^ as a memory 
time during which an clement accumulates kicks before 
attempting a rearrangement. The number of kicks ac- 
cumulated is then proportional to T /Tq. If individual 
kicks are thought of as independent Gaussian perturba- 
tions, and we identify x with the mean-squared size of 
the "cumulative" kick, then x = AT /Tq. The propor- 
tionality constant A would depend, for example, on how 
kicks propagate through the system. For F/Fq — 1, each 
element yields once (on average) within a time interval 
F^ ; A can therefore be viewed as the average number of 
kicks caused by a rearrangement. We leave the analysis 
of such an approach for future work; preliminary investi- 
gations suggest the emergence of interesting features such 
as bistable solutions for the flow curve 17(7). 



B. Attempt frequency To 

Consider now the attempt frequency Fq. It is the only 
source of a characteristic timescale in our model (chosen 
as the time unit above). This excludes a naive proposal 
for the origin of Fq: The attempt frequency cannot be de- 
rived (in some self-consistent way) from the yielding rate 
F, because the model would then no longer contain an in- 
trinsic timescale. This would imply that all dependencies 
on frequency or time are trivial, leading to unphysical re- 
sults (the flow curves (j(^) would simply be a constant, 
as would be the linear moduli G'{uj) and G"{uj)). 

We have so far approximated Fq by a constant value, 
independently of the shear rate 7; this implies that Fq is 
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not caused by the flow directly. One possibility, then, 
is that To arises in fact from true thermal processes, 
i.e., rearrangements of very "fragile" elements with yield 
energies of order k^T. To a first approximation, such 
processes could be accounted for by extending the basic 
equation of motion (^ to 

|p(i?, /, t) = -7|P - r,h e-(^-^'='=)/'=B^ P 

-To e-(^-5fc'')/- p + r(t) piE)dil) (32) 

Here Fth is an attempt rate for true thermal processes, 
which should be a local diffusion rate. In emulsions with 
/xm droplets, typical rates for such diffusive modes could 
be of the order of 1-100 Hz ||ll[. The term on the rhs 
of (^) proportional to Fth corresponds to yield events 
caused directly by thermal fluctuations. Due to the pres- 
ence of interactions between the different elements of 
the material, the effects of such yield events can propa- 
gate through the system and cause other rearrangements. 
These are described by the term proportional to Fq. The 
"attempt frequency" Fq is now no longer an independent 
parameter; instead, it is proportional to the average rate 
of thermal rearrangements, 

The "propagation factor" A again represents the num- 
ber of kicks caused by a thermally induced yield event. 
It has a crucial effect on the behavior of the modified 
model (p2|), as can be seen by considering the equilib- 
rium distribution in the absence of macroscopic strain 
(7(t) = 0). One has PeqiEJ) = Pcq{E)S{l) with 



Poc^iE) = 



Fthe--E/fcBT ^_ pg g-B/x 



p{E) 



When Fq is of the order of Fth or larger, Poq(^) 



exp(£'/x)p(£') as in the original version (g|) of the model. 
From this, the value of Fg can be calculated; for the as- 
sumption Fq ^ Fth to be self-consistent, one then re- 
quires 



Fq _ ^^ JdE p{E)eM~E/kBT) ^ ^ 
Fth JdE p{E) e^E/x) ^ 



(33) 



(here we have neglected a term E/x in the exponent of 
the numerator because k^T <^ x). This condition can be 
given an intuitive interpretation: A must be large enough 
for each thermal yield event to produce at least one new 
element which can yield thermally (i.e., whose yield en- 
ergy E is of order k^T), thus maintaining the popula- 
tion of such fragile elements. For smaller A, one finds 
instead that Fo/Fth ^ exp(— iJ/fceT), which for typical 
barrier energies E = 0{1) (in our units) is unfeasibly 
slow. The above mechanism can therefore give a plau- 
sible rheological time scale only if the average number 
A of rearrangements triggered by one local, thermally 



induced rearrangement is large enough to sustain the 
population of fragile elements, as determined by (|33|). 
The values of A actually required for this are sensitive 
to the small E behavior of p{E). Assuming for example 
p{E) (X E^^^ exp(— i?), one has the condition 



A > [kBTil 



-')]- 



For y = 1, where p{E) stays finite for i5 — > 0, this re- 
quires at least A > 10"'. Such large values appear im- 
plausible unless a single yield event could trigger a whole 
"avalanche" of others; in foams, it has been argued that 
this might be the case pl|. On the other hand, signif- 
icantly smaller values of A would be sufficient if p{E) 
shows a significant bias towards small yield energies E 
(0 « 2/ < 1). The above "thermal trigger" scenario would 
then be more generically plausible. To draw more definite 
conclusions on this point, it would be useful to measure 
p{E) in, for example, a computer simulation of a model 
SGM. 

There are a number of other possible explanations for 
the origin of Fq. These include, for example, noise sources 
internal to the material, such as coarsening in a foam, 
or uncontrolled external noise. Finally, the rheometer 
itself could also be a potential source of noise; this would 
however suggest at least a weak dependence of Fq on the 
shear rate 7. We cannot at present say which of these 
possibilities is most likely, nor rule out other candidates. 
The origin of Fg may not even be universal, but could be 
system specific. 



VII. CONCLUSION 

We have solved exactly the SGR (soft glassy rhcology) 
model of Ref. [|^ for the low frequency shear rheology 
of materials such as foams, emulsions, pastes, slurries, 
etc. The model focuses on the shared features of such 
soft glassy materials (SGMs), namely, structural disor- 
der and metastability. These are built into a generic de- 
scription of the dynamics of mesoscopic elements, with 
interactions represented by a mean-field noise tempera- 
ture x. All rheological properties can be derived from an 
exact constitutive equation. 

In the linear response regime, we found that both the 
storage modulus G' and the loss modulus G" vary with 
frequency as w^"^ for 1 < a; < 2. Near the glass tran- 
sition, they become fiat, in agreement with experimen- 
tal observations on a number of materials. In the glass 
phase, the moduli are predicted to age; this could provide 
an interesting experimental check of the model. 

Far above the glass transition, the steady shear behav- 
ior is Newtonian at small shear rates. Closer to the tran- 
sition (1 < a; < 2), we found power law fluid behavior; in 
the glass phase, there is an additional nonzero yield stress 
(Herschel-Bulkley model) . The last two regimes therefore 
capture important features of experimental data. Above 
the glass transition, the validity of the Cox-Merz rule 
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relating the frequency dependence of the Hnear moduh 
to the shear viscosity can be checked; it breaks down 
in the power law fluid region and fails spectacularly at 
the glass transition. In this regime, stress overshoots in 
shear startup are strongest. We have also calculated the 
distribution of energies dissipated in local yield events. 
At variance with existing simulation data for foams, this 
exhibits a shear-rate dependent crossover between two 
power-law regimes; this discrepancy remains to be re- 
solved. 

We further probed the nonlinear behavior of the model 
by considering large amplitude single and double step 
strains. The nonlinear response cannot in general be 
factorized into strain and time dependent terms, and 
is not well represented by the BKZ equation. Finally, 
we considered measurements of G' and G" in oscillatory 
strain of finite amplitude 7. Near the glass transition, 
G" exhibits a maximum as 7 is increased (strain sweep), 
reproducing qualitative features of recent measurements 
on emulsions and colloidal glasses. The contribution of 
higher harmonics to the stress response is not always suf- 
ficient to determine whether the response is nonlinear, 
emphasizing the need for explicit strain sweeps to get 
reliable data in the linear regime. Otherwise, measure- 
ments at constant strain amplitude can lead to strongly 
enhanced values of the loss modulus G" . Finally, we con- 
sidered the size of the linear regime itself, i.e., the largest 
strain amplitude 7c at which the measured values of G" 
and G" still represent the linear response of the system. 
The SGR model predicts that 7c should be roughly fre- 
quency independent near the glass transition; this point 
should also be amenable to experimental verification. 

In the final section, we speculated on the physical 
origin of the most important parameters of the model, 
namely, the effective temperature x and the attempt fre- 
quency for rearrangements Fq. We argued that x should 
be generically of order unity (in our units). This is be- 
cause it represents the typical energy released in a re- 
arrangement, which is of the same order as the activa- 
tion energy required to cause a rearrangement elsewhere 
in the material. A speculative analogy to marginal dy- 
namics in other glassy systems suggests that x may in 
fact be close to unity in general. This is encouraging, 
because the SGR model reproduces the qualitative fea- 
tures of experimental data best for x « 1, i.e., near the 
glass transition. We mentioned several hypotheses for 
the origin of the attempt frequency Fq, which include 
events triggered by thermal fluctuations or internal and 
external noise sources not explicitly contained within the 
model. 

In future work, we plan to explore in more detail the 
strongly history-dependent behavior of the model in the 
glass phase. Its simplicity should allow this to be done in 
detail, thereby providing the first full theoretical study 
to be made of the generic relationship between aging 
and rheology |4^. Apart from this, the main challenge 
is to incorporate spatial structure and explicit interac- 
tions between elements into the model. This should help 



us understand better the mutual dynamical evolution of 
the attempt rate, the effective noise temperature and the 
structural disorder. In the end, one would hope to de- 
rive a model similar to the present one from such a more 
microscopic description within some well-defined approx- 
imation scheme. 

Acknovifledgements: The author is indebted to 
F. Lequeux, P. Hcbraud and M. E. Gates for their sig- 
nificant contributions to the development and initial in- 
vestigation of the SGR model, as published in [|l^, and 
for helpful comments on the present manuscript. Thanks 
are due also to J. -P. Bouchaud for several seminal sug- 
gestions. Financial support from the Royal Society of 
London through a Dorothy Hodgkin Research Fellowship 
and from the National Science Foundation under Grant 
No. PHY94-07194 is gratefully acknowledged. 



APPENDIX A: DERIVATION OF 
CONSTITUTIVE EQUATION 

The equation of motion (g) of the SGR model can be 
solved by making the time-dependent change of variable 
Z ^ AZ = Z — 7(t). This eliminates the 7 (convective) 
term, converting the equation of motion from a PDF to 
an ODE. Suppressing the E and AZ dependence of P, the 
result reads 



d ( 1 



E-^{Ai+^{t)r 



Pit) 



+r{t)p{E)d{Al + ^{t)) 
This can be integrated to give 



P(t) = P(0) exp 



,^E/, 



z{t,0;Al) 



+ / dt'Tit')piE)6{Al+-fit')) 



'0 

xexp[-e--^/^z(<,f';AOJ (Al) 

with the auxiliary function 

z{t,t';Al)= f dt"exp{[Al + j{t")]y2x} 
Jf 

To simplify matters, we now assume that the initial 
(t = 0) state is completely unstrained, i.e., 7(0) = 
and P(0) ^ Po{E)5{l) = PQiE)S{Al). The stress can 
be calculated by multiplying (Al) by Al and integrating 
over E and Al: 

<7(0 = lit) + (A/)p(,) 

= lit)~ J dt'T(t')^{t') j dEp{E) 

X exp f-e-^/"z(t, t'- -lit'))] (A2) 
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Here the yielding rate r(i) is still undetermined, but it 
can be got from the con diti on of conservation of prob- 
ability: The integral of (Al) over E and A/ has to be 
equal to unity, hence 



1= I dEPo{E)exp -e~^/^z(i,0;0) 
dt'T{t') f dEp{E) 



IQ 

X exp 



-e-^/^z{t,t';-j{t')) 



(A3) 



To write the results (A2,A3) in a more compact form, 
the auxiliary functions defined in (O) and the abbrevia- 
tion (|l|) 

Zit,t')^zit,t';~jit')) 



^dt"exp[[jit")~j{t')f/2x} 



are used. This yields directly eq. ( |10| ) for the yielding 
rate T{t), while for the stress one obtains 

<^{t) - lit) - f dt' Tit') -lit') Gp{Z{t, t')) (A4) 



This can be expressed in the more suggestive form 
by writing the first term on the rhs as 7(i) times the rhs 
of dlT 



APPENDIX B: ASYMPTOTIC BEHAVIOR OF 

GpiZ) 

In this appendix, we derive the asy mpt otic behav- 
ior ( [l3|) of Gp(z). As explained in Sec. Ill, our choice 
of units a^g = 1 implies p{E) = exp[— _E(1 + f{E))] with 
f{E) -^ for E ^ GO. Hence for any 5 > 0, there exists 
M > such that \f{E)\ < 5 ior E > M. Our strategy 
will be to split the defining integral (11) for Gp{z) into 
two parts, for energies above and below the threshold M 
and to bound these separately. Writing 



Gp{z) = / dEp{E)cxj>(-z 
dEp{E)exp(- 



-Ejx 



-E/x 



the first term on the rhs is trivially bounded by zero 
from below and by exp[— zexp(— M/x)] from above. The 
second term, on the other hand, is bracketed by 



rdii;e-(i±^)^/^exp(-. 

JM ^ 



-E/x\ _ 



-2;(1±5) 



dy y 



;(l±5)-lg-a 



(Bl) 



(the plus and minus sign giving the lower and upper 
bound, respectively). Now consider the behavior of 
Gp{z)z^^'^ for some arbitrary small e > 0. Choose 
5 = e/(2x) and a corresponding M\ then from (Bl) 



Gp(z)z''+^ > xz 



c/2 



dy y 



x+e/2-lg-y 



The integral has a finite limit for z — > oo (it is just a 
Gamma function), and so this lower bound tends to in- 
finity in this limit, proving the first part of (|l3|). The 
second part is demonstrated in a similar fashion: With 
the same choice of 5 for a given e, and again using (Bl), 



Gp{z)z'''' < z^-^cxp (-ze-^'/A 
+xz 



'/2 / dy y 

Jo 



,x-e/2-l^-y 



Again, the integral has a finite limit (assuming e is suf- 
ficiently small, i.e., e < 2a;), and both terms on the rhs 
therefore tend to zero for z — > cx), completing the proof 
ofdl). 



APPENDIX C: FLOW CURVES AND YIELD 
STRESS 

Here we derive the small shear rate behavior of the flow 



curves a{j). As shown in Sec. V A , the stress (7(7) = 
/i(7)//o(7) can be expressed in terms of the functions 



Inii) 



dl rGp(Z[l)) 



(CI) 



The scaling of /„ with 7 can be obtained from the asymp- 
totic behavior of Gp{z). From dl3), it follows that for any 
e > 0, we can choose a zq such that 



< Gp{z) < z 



-x+e 



for z > zo 



(C2) 



Now we use zq to decompose the Z-integral in (CI) into 
the parts with I ^ Z07: 



z,'r= / dirGp{z{i)) 

/o 



Replacing Gp{Z{l)) by its minimum and maximum over 
the integration range, I^ is trivially bounded by 



Gp(Z(zo7)) < 



1 



(^07) 



n+l " 



/< <1 



As 7 -^ 0, the Ihs tends to Gp(zo), so we have the 
scaling /< — 0(7"+-^). To bound /^, we use that 
Z{1) > Z/7 > Zq in the relevant integration range, so 
that the bounds ( P2| ) on Gp can be used. Writing Z{1) 
out explicitly, this gives lower and upper bounds for I^ 
of 
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-Tibe 



For X < 71 + 1 (and e sufficiently small), the outer integral 
has a finite limit for 7^0, and so I^ scales as 7^ up 
to sub-power law factors. For larger values of a;, on the 
other hand, this integral diverges as ^"+i-2:±e_ j> t^gn 
scales as 7"'+i (since both the lower and upper bound 
do), i.e., in the same way a s I^. 

As discussed in Sec. VA, the above scaling properties 
of I^ and I^ prove that the flow curve is a power law 
(J r^ 'y^-i (up to sub-power law factors) in the regime 



1 < a; < 2. In the glass phase {x < 1), 
is that of exponential p{E) (eq. (^). 



the simplest case 
The asymptotic 



behavior of Gp{z) 



then translates directly into 



^n ~ 7^ without sub-power law corrections, and this 
gives the Herschel-Bulkley form (gj) of the flow curve. 
The yield stress (pa) is given by the limit of Ii /Iq for 
7 — > 0, while the power law onset of the additional stress 
arises from the small corrections due to Iq . 

For general p(-B), on the other hand, the sub-power 
law factors in Inij) cause a corresponding weak 7 de- 
pendence of cr(7), which dominates the effect of the 
small correction terms 1^(7). The flow curve there- 



fore no longer has the simple Herschel-Bulkley form (24). 
However, in the examples that we tested numerically 
{p{E) - £'"exp(-£') for n = 1, 2, 3), we found that 
this form still provides a good fit to cr{j) over several 
decades of shear rate 7. Both the exponent and yield 
stress of such a fit are then only effective quantities and 
depend on the range of 7 considered; they are therefore 
no longer directly related to x. In the examples that we 
studied, we always found values of the effective exponent 
significantly below unity. 

The slow sub-power law variation of alj) for gen- 
eral p{E) means that there is, for practical purposes, 
always an effective yield stress (whose actual value de- 
pends weakly on the lowest accessible shear rate 7) . Nev- 
ertheless, one may wonder what the "true" yield stress 
CTy = (7(7 —>■ 0) would be. The above line of argument 
does not answer this question; it does not even exclude 
the possibility of Cy being zero. We have examined this 
issue for several different sub-power law corrections to the 
asymptotic behavior of Gp, such as Gp{z)z^ ~ (Inz)™, 
or ~ exp[(lnz)"] with \n\ < 1. The yield stress is al- 
ways nonzero, and in fact turns out to be the same as for 
exponential p{E). We suspect that this may be true in 
general, but have not found a proof. 



APPENDIX D: NUMERICAL DETERMINATION 

OF G*((.<j,7) 

In this appendix, we outline the numerical scheme 
that we used to obtain the nonlinear dynamic modulus 
G*{ui,j) and the residual r defin ed in ( po|) and (|3^), re- 
spectively. As explained in Sec. VD, we are interested 



in the steady state stress response in the ergodic regime 
X > 1. We can then safely send the initial time to —00 in 
the CE (^,^. The equations that need to be solved can 
be simplified further by using the fact that in the steady 
state, the yielding rate T(t) must have the same period- 
icity as the applied strain j{t). Denoting the oscillation 
period by T = 2tt/u!, the task is then to solve 



dt' T{t')H{t,t') 



t-T 



for r(i) and then to evaluate the stress from 



a{t) = 7(0 



t-T 



dt' -f{t')T{t')H{t,t') 



(Dl) 



(D2) 



Here the periodicity of the problem has been absorbed 
into the definition of 



7I(i,t') = ^Gp(Z(t,i'-nT)) 

= J2 Gp{Z{t, t') + nZ{t' + T, t')) 



n=0 



where the second equality follows again from the period- 
icity of the strain 7(t) = 7C0S wt. The numerical solution 
of the integral equation (Dl) is simplified by subtracting 
from the kernel H{t, t') a part that depends on t' only: 



H{t, t') = H{t, t') - H{t' + T, t') = 



p-nzi _ p-^z-i 



1 



,-OZ2 



where we have abbreviated VL = exp(— i?/x), Z\ = 
Z{t,t'), Z2 = Z{t' + T,t'). The modified kernel H{t,t') 
has the convenient properties H{t' ,t') = 1, H{t'+T,t') = 
and is also simpler to evaluate numerically than H{t, t'). 
The yielding rate can easily be calculated from H instead 
of H: Defining a modified yielding rate V{t) as the solu- 
tion of 



1 = 



t-T 



dt' f{t')H{t,t') 



(D3) 



the actual yielding rate is recovered by dividing by the 
constant factor 



dt' T{t')H{t' + T,t') 



However, even the solution of (p3) is still nontrivial, es- 
pecially in the low frequency regime T ^ I that we are 
most interested in. This is because H "inherits" from 
Gp an initial "fast" decay as t — t' increases from zero, 
followed by a much slower power-law decay (which in 
turns gives way to a rapid final decay as soon as strain- 
induced yielding becomes important). This separation of 
0(1) and 0{T) timescales rules out traditional solution 
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metho ds s uch as Chebyshev approximation. Instead, we 
solve (D3) by Fourier transform: Writing 



n— — oo 

eq. ( p3| ) is transformed into the matrix equation 

oo 
m— — oo 

with coefficients 



(D4) 



Hr, 



dt 
T 



e -^C"-™)'^* / dT e-™'"^iJ(<, t - r) 



Once (D4) is solved and the rescaling from F to F is 
carried out, the stress is obtained as 



a{t) 



Y.(Tnt 



'^n — T^i^n,-! + <5n,l) ~ ^ /_^ ^7n{Hn^,n+l + i?„,m-l) 



Its Fourier components determine the nonlinear dynamic 
modulus and squared residual as 



G*iuj,-f)^2ai r^ 



1 



^1 



EOO I 

fc=0 F2fe4 



The result for r^ has been simplified using the fact that 
(T_„ = (7* (because a{t) is real) and that cr„ = for even 
n (because a{t) -^ — o'(i) for 7 -^ —7, which corresponds 
tot^t + T/2). 

To solve the main equation ( p4| ), we truncate the ma- 
trix equation at successively higher orders until the cal- 
culated values of G"(ci;,7), G"{uj,j) and r are stable to 
within 1%. The Fourier components Hmn are calculated 
from a spline interpolant approximation to H(t, t') in or- 
der to save expensive function evaluations. 
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